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MONTE  CARLO  METHODS  FOR  NEUTRON  FLUX  CALCULATIONS  IN 
A  PRESSURIZED  LIGHT  WATER  POWER  REACTOR  USING  MORSE -CG 

I.  INTRODUCTION 


The  fast  neutron  flux  (En>1.0  MeV)  and  neutron  energy  spectra 
in  the  core  midplane  of  a  typical  Oconee  class  pressurized  light  water 
power  reactor  have  been  calculated  using  the  Monte  Carlo  computer  code 
MORSE1'  with  combinatorial  geometry. 


The  purpose  of  the  calculation  was  three  fold:  (1)  to  compare 
fast  neutron  flux  in  the  core  midplane  of  a  nuclear  reactor  using  a 

one-dimensional  slab  model  as  calculated  with  the  code  MORSE  with  that 

2  3 

calculated  using  the  discrete  ordinates  code  ANISN  ;  (2)  to 

compare  fast  neutron  flux  results  in  slab  and  cylindrical  models;  and 

(3)  to  investigate  the  effect  on  neutron  flux  results  of  using 

different  flux  estimators. 


Note:  Manuscript  submitted  September  20,  1979. 
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The  comparison  between  MORSE  and  ANISN  using  the  slab  model  was 
investigated  both  to  use  as  a  benchmark  calculation  and  to  check  the 
accuracy  of  ANISN  against  MORSE  in  a  reactor  geometry  situation. 

ANISN  solves  the  Boltzman  transport  equation  in  a  discretized  space 
with  the  method  of  finite  differences,  while  MORSE  solves  the  Boltzman 
transport  equation  in  a  continuous  space  using  a  Monte  Carlo  method. 
Although  these  two  methods  have  been  compared  previously  for  simple 
homogenous  slabs,  the  accuracy  of  ANISN  when  applied  to  amultiple-slab 
reactor  geometry  has  not  been  checked  in  detail.  Comparison  between 
the  slab  and  cylindrical  models  was  also  made. 

In  Monte  Carlo  type  calculations,  it  is  generally  of  some  concern 
as  to  which  estimator  to  use.  An  estimator  is  generally  selected  that 
will  minimize  the  number  of  histories  needed  to  attain  a  given 
statistical  accuracy.  It  is,  therefore,  useful  to  compare  fast 
neutron  flux  results  in  the  slab  model  using  various  types  of 
estimators.  For  the  slab  and  cylindrical  geometries,  the  only  flux 
estimators  we  seriously  considered  were  the  boundary  crossing 
estimator  and  the  next  event  uncollided  flux  estimator.  We  have 
therefore  compared  results  for  these  two  estimators,  but  also  used  the 
collision  density  and  track  length  per  unit  volume  estimators4  as 
checks. 
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The  slab  geometry  and  neutron  source  parameters  used  to  describe 
the  reactor  are  discussed  in  detail  in  Section  II.  Section  III  is  a 
discussion  of  the  biasing  of  the  sampling  of  the  neutron  source  energy 
and  spatial  distributions  in  the  reactor  core.  In  addition 
path-length  stretching,  energy  biasing  at  the  collision  site,  Russian 
roulette,  splitting,  and  angle  biasing  are  discussed  there.  Various 
neutron  flux  estimators  are  discussed  in  Section  IV  and  the 
cylindrical  modeling  of  the  reactor  in  Section  V.  A  discussion  of  the 
results  of  the  fast  neutron  flux  and  neutron  energy  spectra  calculated 
by  MORSE  and  ANISN  are  given  in  Section  VI. 

II.  MORSE  MONTE  CARLO  CALCULATION  IN  THE  SLAB  MODEL 

A.  Composition  and  dimensions  of  reactor  components 

The  characteristics  of  the  reactor  in  this  calculation  were 

2 

typical  of  Babcock  and  Wilcox  Oconee  class  power  reactors.  Neutron 

2 

flux  and  spectrum  calculations  using  ANISN  have  been  performed  by 
3abcock  and  Wilcox  using  these  reactor  parameters  and  the  cross 
section  library  CASK.^ 

In  the  slab  model  each  component  or  region  of  the  reactor  is 
treated  as  a  infinite  slab  with  the  normal  to  the  surface  in  the  Z 
direction  as  shown  in  Fig.  1.  For  the  MORSE  calculation  the  size  of 
the  slab  in  the  X  and  Y  direction  was  actually  taken  to  be  20,000  cm, 


3 


which  is  infinite  for  all  practical  purposes.  The  exact  dimensions  of 
each  slab  is  given  in  Table  1.  The  elemental  composition,  and  atomic 
densities  in  each  region  are  given  in  Table  2.  The  atomic  densities 
in  the  reactor  core  were  determined  assuming  the  core  is  a  homogeneous 
mixture. 

B.  Core  power  distribution  function 

The  relative  core  power  distribution  was  obtained  from 

criticality  calculations  by  Babcock  and  Wilcox  using  the  computer 

6  7 

codes  ’  PDQ-5  and  Harmony.  In  the  slab  model  the  source  strength 
is  proportional  to  the  product  of  a  relative  power  distribution 
function  P(Z),  which  is  invariant  under  translations  in  the  X  or  Y 
direction  and  a  neutron  energy  distribution  function  x (E) .  We  assume 
that  source  neutrons  are  emitted  isotropically.  The  function  P(Z)  is 
given  as  a  discrete  set  of  average  values  on  a  set  of  intervals  in  the 
core  region  and  are  tabulated  in  Table  3  (together  with  the  relative 
power  density).  The  energy  distribution  function  x(E)  is  similarly 
defined  on  a  discrete  set  of  energy  intervals  or  groups.  The  cross 

(4  \ 

section  library  used  in  the  calculation  is  CASK,'  '  a  40  group 
coupled  neutron  and  gamma-ray  cross-section  data  set  (22  neutron  and 
18  gaima-ray  groups).  We  only  made  use  of  the  14  nighest  energy 
neutron  groups  since  we  are  interested  in  fast  flux  only.  The 
function  x  is  tabulated  in  Table  4.  The  neutron  source  strength 
function  can  then  be  written  as 


Table  1 


Region  dimensions  of  components  of  the  slab  geometry  reactor 


REGION 

REACTOR 

COMPONENT 

OUTER 

BOUNDARY 

(cm) 

REGION 

THICKNESS 

(cm) 

1 

CORE 

163.58 

163.58 

2 

COOLANT 

163.79 

0.21 

3 

LINER 

165.70 

1.91 

4 

COOLANT 

179.07 

13.37 

5 

CORE  BARREL 

184.15 

5.08 

6 

COOLANT 

186.69 

2.54 

7 

SUPPORT 

CYLINDER 

191.77 

5.08 

8 

COOLANT 

217.17 

25.40 

9 

PRESSURE 

VESSEL 

238.44 

21.27 

10 

CAVITY 

350.52 

116.08 

11 

PRIMARY  SHIELD 

502.92 

152.40 
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Table  2 

Material  composition  of  the  reactor  components 


REGION 


MATERIAL 


ATOMIC  DENSITY 
atoms/ cm3 


7 


10 


n 


I 


PRESSURE  VESSEL  STEEL  (A5338) 


CARBON 

8.6704E 

+ 

20 

ALUMINUM 

7.0177E 

+ 

19 

SILICON 

4.2029E 

+ 

20 

CHROMIUM 

1.2746E 

+ 

20 

MANGANESE 

1.1201E 

+ 

21 

IRON 

3.1974E 

+ 

22 

NICKEL 

4.8377E 

+ 

20 

MOLYBDENUM 

2.7137E 

+ 

20 

CAVITY 


AIR 

( 150°F 


15PSI) 

NITROGEN 

OXYGEN 


3.413E  +  19 
9.156E  +  18 


PRIMARY  SHIELD  ORDINARY  CONCRETE 

HYDROGEN 

CARBON 

OXYGEN 

SODIUM 

MAGNESIUM 

ALUMINUM 

SILICON 

POTASSIUM 

CALCIUM 

IRON 


8.6089E  +  21 
1.1423E  +  20 
4.3289E  +  22 
9.6396E  +  20 
1.2392E  +  20 
1.741 OE  +  21 
1.6618E  +  22 
4.6060E  +  20 
1.5028E  +  21 
3.4503E  +  20 
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(1) 


S(Z,E)  =  1.323  x  105  P(Z)  (E)  n/sec 

The  proportionality  constant  was  chosen  to  express  the  strength 
function  for  this  reactor  in  units  of  neutrons  per  second.  It  also 
includes  a  correction  factor  of  1.55  (the  product  of  azimuthal  and 
axial  peaking  factors)  to  take  into  account  the  fact  that  the  maximum 
flux  occurs  out  of  the  mid-core  plane  since  the  core  and  its  power 
distribution  are  really  not  slab  or  cyli ndr ical ly  symmetric. 

Because  of  the  symmetry  of  the  core  power  distribution  and  the  reactor 
medium,  one  need  not  sample  the  entire  volume  source,  but  only  the 
line  source  distribution  along  X=Y=0.  The  thickness  of  the  core  slab 
was  taken  equal  to  the  core  radius;  a  specular  reflection  boundary 
condition  at  the  plane  Z=0  was  used  to  more  closely  correspond  to  the 
near  cylindrical  geometry  of  the  core. 

III.  3IASING 

A.  Spatial  biasing  of  the  source 

Since  the  core  has  a  thickness  of  163.58  cm  and  we  wish  to 
know  the  flux  out  to  and  through  the  pressure  vessel  wall  which  is 
many  mean  free  paths  away,  it  is  necessary  to  bias  the  source  spatial 
distribution  so  that  most  of  the  neutrons  are  selected  from  near  the 
surface.  Otherwise  few  of  the  source  neutrons  selected  will  get  out 


Table  3 

Core  power  distribution  for  slab  geometry 


INTERVAL  INTERVAL 

MIDPOINT  WIDTH 

(cm)  (cm) 


RELATIVE 
POWER  DENSITY 


( cm" 1 ) 


SLAB 

RELATIVE 

POWER 

P(z) 


8.3333 

25.000 

41.667 
58.333 
75.000 

91.667 

101.36 
104.09 
106.32 
109.55 
112.27 
115.00 
117.73 

120.45 
123.18 
125.91 
128.64 

130.54 

131.62 

132.71 

133.79 

134.87 

135.96 
137.04 

138.12 

139.21 

140.29 

141.37 

142.46 

143.54 

144.62 

145.71 

146.79 

147.87 

148.96 
150.04 

151.12 

152.21 

153.29 


16.6667 

II 

II 

II 

II 

II 

2.72727 

II 

II 

II 

II 

II 

II 

II 

II 

II 

II 

1.08323 

It 

II 

II 

ll 

II 

II 

II 

II 

II 

II 

II 

II 

II 

II 

II 

II 

II 

II 

II 

1.08323 

II 


0.94226 

.97196 

.99001 

.97394 

.94643 

.95221 

.96924 

.96529 

.95953 

.95569 

.96169 

.97221 

.98164 

1.0095 

1.0672 

1.1574 

1.2350 

1.2581 

1.2913 

1.3170 

1.3370 

1.3407 

1.3291 

1.3091 

1.2870 

1.2665 

1.2495 

1.2362 

1.2264 

1.2190 

1.2120 

1.2034 

1.9111 

1.1732 

1.1516 

1.1303 

1.1047 

1.0714 

1.0304 


15.705 
16. 198 
16.500 
16.232 
15.773 
15.370 
2.6434 
2.6326 
2.6169 
2.6064 
2.6234 
2.6542 
2.6772 
2.7532 
2.9105 
3.1565 
3.3682 
1.3628 
1.3988 
1.4266 
1.4483 
1.4523 
1.4397 
1.4180 
1.4049 
1.3719 
1.3535 
1.3391 
1.3284 
1.3205 
1.3129 
1.3036 
1.2902 
1.2708 
1.2474 
1.2244 
1.1966 
1.1606 
1.1162 
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154.37 

.97903 

1.0605 

155.46 

.92429 

1.0012 

156.54 

.86173 

.9335 

157.62 

.80526 

.87228 

158.71 

.74459 

.80656 

159.79 

. 70086 

.75919 

160.87 

.66114 

.71617 

161.96 

.62466 

.67665 

163.04 

.58864 

.63763 
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Table  4 


Source  neutron  energy  spectrum 


ENERGY  GROUP 
9 

NEUTRON  ENERGY 

NO.  UPPER  EDGE  GROUP  ENERGY 
(MeV) 

NORMALIZED 

RELATIVE  PROBABILITY 

(g) 

1 

14.92 

0. 16000E-03 

2 

12.20 

0.90000E-03 

3 

10.00 

0.35000E-02 

4 

8.180 

0. 13970E-01 

5 

6.360 

0.34730E-01 

6 

4.960 

0. 35220E-01 

7 

4.060 

0. 10778E+00 

8 

3.010 

0.89410E-01 

9 

2.460 

0.23300E-01 

10 

2.350 

0. 12091E+00 

n 

1.830 

0.2 1913E+00 

12 

1.110 

0. 19937E+00 

13 

.5500 

0. 13605E+00 

14 

.1110 

0.15570E-01 

12 


and  an  inordinate  number  of  histories  will  be  needed  to  obtain  a 

reasonable  variance  in  the  flux  at  the  pressure  vessel.  The  original 

neutron  spatial  distribution  ?rZ)  was  biased  therefore  by  the  factor 

expf (Z  -Z)/ X  ) ,  where  Z  is  the  core  radius  and  \  is  an  averaae 
o  o 

neutron  mean  free  path  in  the  core,  so  that  the  frequency  of  selection 
of  the  neutron  starting  position  in  the  core  was  roughly  oroportional 
to  the  probability  of  the  neutron  escaping  from  the  core.  The  weight 
of  the  selected  neutron  is  then  modified  so  that  for  each  interval  the 
product  of  starting  weight,  W.,  and  the  probability  of  Dicking  a 
neutron  from  that  interval  is  invariant.  That  is,  if  Z.  is  the 
midpoint  of  the  ith  interval 

WiP(Zi)  =  Cwjp'tf.)  =  cw'p(z.)  exp  -fZ0  -  Zi)/X] 
where  C=  SPfZ^exp  (Z1-  -  ZQ ) /X ] 

» 

Therefore,  the  modified  weight  of  the  neutron,  A.,  is 

w’.  =  (Wi/C)  exp[  +  (Z0-Zi )/X]  (4) 

In  the  actual  calculation  we  took  x  =  7  cm.  The  MORSE  subroutine 
SOURCE  was  'written  to  incorporate  this  biasing  feature  and  is  listed 
in  Appendix  I. 


(21 

(3) 
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3.  Source  energy  biasing 


The  source  neutron  energy  distribution  function  x(E)  as 
shown  in  Table  4  Deaks  in  the  1-2  MeV  range  and  the  neutrons  in  the 
10-15  Me'/  range  are  100-1000  times  less  orobable.  It  was  necessary  to 
emoloy  energy  biasing  in  the  source,  since  we  are  interested  in  the 
fast  neutron  flux  out  to  the  pressure  vessel  wall  and  contributions  to 
the  flux  will  be  greater  for  those  neutrons  starting  out  with  the 
higher  energies.  The  lower  energy  neutrons  will  undergo  so  many 
collisions  that  their  energy  will  be  well  below  1  MeV  before  they 
reach  the  wall  and  their  contributions  to  the  fast  flux  will  be 
negligible.  An  energy  biasing  function  3(E)  was  chosen  such  that  the 
frequency  of  selection  of  neutron  energies  was  apDroximately  uniform 
over  the  energy  region  sampled.  That  is  all  energy  groups  are  equally 
likely  to  be  chosen.  In  this  case  the  modified  energy  distribution, 
x'(E^)  is  obtained  by  multiplying  the  corresponding  x(E.)  by  the 
biasing  factor  b  /8 ( E ,  where  3/E)  is  tabulated  in  Table  5,  and  the 
constant,  b,  is  such  that  the  new  distribution  x'  is  normalized. 

C.  Path-length  stretching 

Another  form  of  biasing  that  we  employed  in  the  calculation 
of  the  neutron  flux  is  known  as  path-length  stretching.  Instead  of 
using  the  physical  mean  free  path,  X,  between  collisions,  the  mean 
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Table  5 


Energy  viasing  factors  for  source  neutrons  and  at  collision  sites 


GROUP  No. 

UPPER  EDGE 
ENERGY  GROUP 
MeV 

SOURCE  NEUTRON 
ENERGY  BIASING 
FACTOR 

BfE) 

COLLISION 
ENERGY  3IASING 
FACTOR 

1 

14.92 

0.0495 

1.28  E+02 

2 

12.20 

0.1052 

6.40  E+01 

3 

10.00 

0.1702 

3.20  E+01 

4 

3.180 

0.2481 

1.60  E+01 

5 

6.360 

0.3125 

8.00  E+00 

6 

4.960 

0.3780 

4.00  E+00 

7 

4.060 

0.4447 

2.00  E+00 

3 

3.010 

0.5001 

1.75  E+00 

9 

2.460 

0.5722 

1.50  E+00 

10 

2.350 

0.6471 

1.25  E+00 

n 

1.830 

0.7827 

1.00  E+00 

12 

1.110 

0.9061 

7.50  E-01 

13 

0.550 

0.9904 

6.00  E-02 

14 

0.110 

1.000 

3.00  E-02 

L5 


free  path  in  a  direction  of  interest  can  be  effectively  increased  by 
factor  denoted  "BIAS"  in  MORSE,  so  that,  e.g.,  the  neutron  travels 
further  in  the  direction  of  the  pressure  vessel  before  encountering  a 
collision. 


The  factor  BIAS  is  defined  in  terms  of  two  other  Darameters, 


where  PATH  is  a  measure  of  how  much  stretching  is  to  be  applied  and 
lies  between  0  and  1  and  DIREC  is  taken  as  the  cosine  of  the  angle 
between  the  flight  direction  and  the  direction  in  which  one  wants  to 
encourage  the  neutrons  to  go.  This  technique  will  improve  counting 
statistics  in  regions  in  the  direction  of  the  pressure  vessel  at  the 
expense  of  a  loss  in  statistics  in  the  local  region  where  the  biasing 
is  applied.  When  properly  apolied  this  tradeoff  is  beneficial.  The 
results  tabulated  herein  were  calculated  with  path  stretching 
parameters  (PATH)  shown  in  Table  5. 

0.  Energy  biasing  at  collision  sites 

As  a  neutron  enters  a  collision  its  outgoing  energy  and 
direction  after  collision  are  determined  by  first  sampling  the  group 
to  group  transfer  matrices  to  obtain  the  new  energy  and  then  the 


Table  6 


Values  of  the  path  length  stretching  parameter  PATH 


REGION 

REACTOR 

COMPONENT 

VALUE  OF 
"PATH" 

1 

Core 

0.75 

2 

Coolant 

0.0 

3 

L  i  ner 

0.0 

4 

Coolant 

0.5 

5 

Core  Barrel 

0.0 

6 

Coolant 

0.0 

7 

Support  Cylinder 

0.0 

8 

Coolant 

0.5 

9 

Pressure  Vessel 

0.5 

10 

Cavity 

0.0 

11 

Primary  Shield 

0.6 
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angular  distribution  that  is  kinematically  consistent  with  this 
neutron  energy.  Low  energy  neutrons  near  the  core  are  less  likely  to 
contribute  to  the  fast  flux  at  the  vessel  wall,  and  so  one  would  like 
to  bias  the  outgoing  energy  distribution  toward  higher  energies. 

In  MORSE  one  can  also  bias  the  sampling  of  the  group  to 
group  energy  transfer  probability  at  collision  sites  so  that  by 
employing  these  multiplicative  bias  factors,  on  the  average  higher 
energy  neutrons  are  selected  at  each  collision.  The  set  of  bias 
factors  used  in  this  calculation  is  shown  in  Table  7. 

The  comparison  in  slab  and  cylindrical  geometry  between 
MORSE  and  ANISN  in  this  report  does  not  include  collision  biasing. 

E.  Russian  roulette  and  splitting 

There  are  two  other  options  available  in  MORSE  for 
decreasing  variance  and  increasing  efficiency.  When  the  weight  of  a 
neutron  becomes  so  small  that  it  is  inefficient  to  follow  its  history, 
one  can  "play"  Russian  Roulette.  With  a  certain  probability,  the 
particle  is  either  "killed"  or  its  weight  is  increased  so  that  it  pays 
to  follow  it  again.  When  the  weight  of  a  particle,  on  the  other  hand, 
is  too  large  and  its  contribution,  or  lack  of  it,  may  cause  a  large 
fluctuation  in  the  final  result,  one  may  split  that  particle  into  a 
number  of  particles  each  having  a  smaller  weight.  We  found  neither  of 
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Table  7 

Comparison  of  fast  neutron  flux  using  MORSE  and  ANISN  in  slab  geometry 


RADIUS 
( cm) 

ANISN 

(Flux) 

E  >  1 .0  MeV 

MORSE (SLAB) 
(Flux) 

E  >  1 .0  MeV 

Standard 

Deviation 

% 

M0RSE/ANISN 

163.58 

4.60E  +13 

4.58E  +13 

+2 

1.00 

163.79 

4.42E  +13 

4.39E  +13 

+2 

0.99 

165.70 

2.64E  +13 

2.62E  +13 

+2 

0.99 

169.04 

1.35E  +13 

1.43E  +13 

+3 

1.06 

172.39 

8.27E  +12 

8.61E  +12 

+3 

1.04 

175.72 

5.49E  +12 

5.70E  +12 

+4 

1.04 

179.07 

4.29E  +12 

4.52E  +12 

+4 

1.05 

181.61 

3.25E  +12 

3.12E  +12 

+4 

0.96 

134. 15 

1.98E  +12 

1.84E  +12 

+4 

0.93 

186.69 

1.43E  +12 

1.45E  +12 

+5 

1.01 

191.77 

5.93E  +11 

6.21E  +11 

+  13 

1.05 

200.24 

1.48E  +11 

1.64E  +11 

+  11 

1.11 

208.71 

5.36E  +10 

5.98E  +10 

+34 

1.12 

217.17 

2.80E  +10 

2.03E  +10 

+  15 

0.73 

222.49 

1.68E  +10 

1.09E  +10 

+20 

0.65 

227.81 

8.36E  +9 

1.09E  +10 

+42 

1.13 

233.12 

4.00E  +9 

3.97E  +9 

+25- 

.99 

238.44 

1.73E  +9 

1.66E  +9 

+38 

0.96 
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these  options  to  be  very  useful  in  our  calculation  of  fast  flux  and  so 
did  not  use  them. 

IV.  NEUTRON  FLUX  ESTIMATORS 


A.  3oundary  crossing  or  surface  crossing  estimators  (BDRYX) 


This  method  is  particularly  useful  for  estimating  the  flux  in 
a  one  dimensional  geometry.  We  can  obtain  the  contribution  to  the 
average  flux  on  each  boundary  or  surface  crossed  by  the  ith  neutron  by 
scoring  the  weight  (when  crossing)  of  that  neutron,  W^,  divided  by 
the  absolute  value  of  the  cosine  of  the  angle  between  the  neutron 
track  and  the  normal  to  the  surface: 


♦  (Z) 


.£  A_ 

i  |cos  9  .j 


(6) 


Because  of  the  symmetry  of  the  medium  and  the  source 
distribution,  the  flux  is  uniform  on  any  symmetry  surface.  Hence  the 
flux  at  any  point  on  that  surface  is  equal  to  the  average  flux  over 
the  surface. 


A  difficulty  with  the  boundary  crossing  estimator  is  that  the 
variance  associated  with  it  is  unbounded.  It  is  clear  that  for 
grazing  angles  9  ~  w/2,  1/cos  9  is  very  large,  making  a  very  large 
contribution  to  the  sum  and  the  variance.  To  avoid  these  large 
fluctuations,  we  limit  the  magnitude  of  cos  9  ,  by  setting  cos  9  to 
0.005  for  all  angles  such  that  cos  9  <  .01. 
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This  estimator  is  implemented  in  the  subroutine  BDRYX  and  is 
listed  in  Appendix  II. 

The  subroutine  BDRYX  is  called  whenever  a  neutron  track 
crosses  the  boundary  of  a  geometry  medium.  Even  if  the  problem  has 
only  one  real  medium,  different  geometry  media  can  be  artificially 
defined  so  that  scoring  surfaces  can  be  established  wherever  they  are 
desired.  These  geometry  media  are  defined  in  the  input  to  MORSE  and 
obtained  from  subroutine  GTMED,  which  is  listed  in  Appendix  III. 

3.  Next-event  uncollided  estimator  (UNC) 


In  the  surface  crossing  estimator  previously  discussed 
neutrons  had  to  actually  cross  the  surface  or  boundary  in  order  to  be 
scored.  We  can  also  use  a  so  called  expected  value  estimator  which 
relates  the  emergent  particle  density  at  a  collision  site  to  the 
flux  $(Z).  Thus  a  contribution  to  the  flux  is  made  at  every 


collision.  The  flux  as  calculated  by  this  estimator  is  given  by 

+  U(Z)  (7) 


^  W  .  exp  (-  E.R., 

$  (Z)  =  Z  — 1 - ^ 

i 


COS  9  • 


where  W..  is  the  statistical  weight  after  collision,  z t ( E ^ )  is 
the  macroscopic  total  neutron  scattering  cross  section  for  the 
emergent  particle  of  energy  E.,  R..  is  the  distance,  in  the 
direction  of  the  velocity  of  the  emergent  particle,  between  the  ith 
collision  site  and  the  plane.  This  flux  estimate  is  implemented  in 


21 


subroutine  RELCOl.  3ecau$e  this  subroutine  is  only  called  at 
collision  sites,  a  separate  analytic  contribution,  UfZ),  from  the 
source  site  must  be  made — the  uncollided  flux  contribution. 

C.  Tracklength  per  unit  volume  estimator  (TLPUV) 

The  neutron  weighted  total  path  length  within  some  control 

volume,  divided  by  that  volume  is  the  average  flux  throughout  the 

.  4 

volume: 

$  =  —  £w  L  (8) 

4  V  i  i  i 

A  large  control  volume  can  be  chosen  to  improve  statistics,  but  a 
small  control  volume  will  give  a  more  accurate  value  of  the  flux  at  a 
point.  We  have  also  included  this  estimator  for  comparison  with  the 
others.  In  MORSE  the  estimate  can  be  obtained  from  subroutine 
ENDR'JN.  However  no  variance  on  the  estimate  is  given  there. 

If  the  neutron  suffers  many  collisions  within  the  control 
volume,  this  method  becomes  rather  inefficient  and  estimators  based 
on  the  density  of  collisions  should  be  considered. 
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D.  Collision  density  concept  of  flux  (COLDEN) 


If  within  some  control  volume,  AV,  the  weighted  total  number 

of  Monte  Carlo  collisions  of  neutrons  in  energy  group  g-  is  given  by 

4 

n(g^),  then  the  average  flux  in  the  control  volume  is  given  by 

"  l9;>  .  (9) 

-  i  Et(g.)AV 

where  2  is  the  total  macroscopic  scattering  cross  section.  This 
method  works  well  when  there  are  many  scattering  collisions  in  the 
control  volume.  Again,  for  comparison  purposes,  a  flux  estimate  based 
on  this  model  is  also  given.  It  also  can  be  obtained  from  the 
subroutine  ENDRUN. 

V.  VARIANCE  REDUCTION  VIA  PATHLENGHT  STRETCHING  AND  ENERGY  BIASING  AT 
COLLISION  SITES  (Collision  Biasing) 

In  the  simple  slab  geometry  we  have  computed  the  fast  flux  as  a 
function  of  radial  distance  from  the  reactor  core  for  900  neutron 
histories  to  determine  how  path  length  stretching  and  energy  biasing 
at  collisions  affect  the  variance  for  the  boundary  crossing  and 
next-event  uncollided  estimators.  These  were  test  cases  performed  to 
look  for  qualitative  trends  and  clearly  were  not  made  to  reduce  the 
variance  to  the  lowest  possible  value. 
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The  calculated  fractional  standard  deviation  for  the  UNC 
estimator  is  plotted  versus  radial  distance  from  the  core  in  Fig.  2 
for  the  path  length  stretching  parameters  shown  in  Table  6.  On  the 
average,  there  is  a  definite  variance  reduction  for  distances  beyond 
180  cm  for  the  path  length  stretched  calculations.  For  distances  less 
than  180  cm  there  are  no  observable  differences  . 

A  similar  comparison  is  made  for  the  UNC  estimator  with  and 
without  collision  biasing.  The  collision  biasing  factors  were 
optimized  empirically  and  are  tabulated  in  Table  6.  A  plot  of  these 
results  is  shown  in  Fig.  3.  Again  there  is  a  significant  variance 
reduction  for  the  fast  flux  beyond  distances  of  190  cm.  For  distances 
less  than  190  cm  we  were  able  to  see  no  significant  difference  in  the 
vari ance. 

In  Fig.  4  we  have  plotted  the  results  when  both  collision  biasing 
and  path  length  stretching  were  used  and  when  neither  was  used. 

Again,  beyond  190  cm,  there  is  a  significant  reduction  in  variance  for 
the  biased  case.  We  conclude  that  utilization  of  collision  biasing 
and  path  length  stretching  will  yield  more  accurate  flux  values  near 
the  pressure  vessel  wall  for  a  given  number  of  neutron  histories. 
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FRACTIONAL  STANDARD  DEVIATION  IN  FAST  FLU 


Fig.  2  -  Fractional  standard  deviation  in  the  midcore  plane  fast 
flux  with  and  without  stretching 


FRACTIONAL  STANDARD  DEVIATION  IN  FAST  FLUX 


FRACTIONAL  STANDARD  DEVIATION  IN  FAST  FLUX 


VI.  COMPARISON  OF  SORYX  WITH  (JNC 


We  have  performed  a  similar  analysis  for  the  boundary  crossing 
estimator  (BDRYX).  The  plots  of  fractional  standard  deviation  versus 
radial  distance  are  shown  in  Figs.  5  and  6  for  collision  biasing  and 
stretching,  respectively.  The  results  in  neither  case  are  clear 
because  an  insufficient  number  of  neutron  histories  was  taken. 

A  comparison  of  the  fractional  standard  deviation  versus  radial 
distance  for  the  estimators  BDRYX  and  UNC  with  collision  biasing  and 
path  length  stretching  included  for  both  cases  is  shown  in  Fig.  7. 
There  seems  to  be  no  significant  difference  anywhere  in  the  outside 
the  core  region.  For  this  simple  slab  geometry  UNC  does  have  the 
advantage  in  that  it  required  about  13  secs  of  computer  CPU  time  while 
the  BDRYX  estimator  required  21  secs  for  the  same  900  histories. 
However,  for  more  complicated  geometries,  this  advantage  may  not 
pers ist. 

VII.  COMPARISON  OF  THE  MONTE  CARLO  SLAB  MODEL  FAST  FLUX  RESULTS  WITH 
ANISN 

The  final  slab  geometry  Monte  Carlo  calculations  were  performed 
using  200,000  neutron  histories  and  the  boundary  crossing  estimator. 

At  the  same  time  the  fast  flux  was  also  estimated  by  track  length  per 
unit  volume  and  collision  density  estimators.  These  last  two 
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D  DEVIATION 


Fractional  standard  deviation  in  the  midcore  plane  fast 
flux  for  BDRYX  estimator  with  collision  biasing 


FRACTIONAL  STANDARD  DEVIATION  IN  FAST  FLUX 


f 


z  (  c  m  ) 

Fig.  6  -  Fractional  standard  deviation  in  the  midcore  plane  fast 
flux  for  BDRYX  estimator  with  stretching 
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FRACTIONAL  STANDARD  DEVIATION  IN  FAST  FLUX 


Fig.  7  -  Fractional  standard  deviation  in  the  midcore  plane  fast  flux 
for  BDRYX  and  UNC  with  collision  biasing  and  stretching 
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estimators  were  useful  as  a  check  on  BORYX  to  be  sure  no  gross  errors 
were  built  into  8DRYX.  Since  these  two  estimators  averaged  flux  over 
a  volume  between  two  regions  instead  of  giving  the  flux  at  a  boundary 
of  a  region,  an  exact  comparison  to  BDRYX  cannot  be  expected.  Source 
energy  and  spatial  biasing  were  used.  The  path  length  stretching 
parameters  were  empirically  optimized  and  the  final  set  is  tabulated 
in  Table  1.  No  energy  biasing  at  the  collision  sites  was  employed.  A 
graph  of  the  relative  flux  estimated  by  the  collision  density 
(COLDEN),  tracklength  (TLPUV)  and  boundary  crossing  estimators  (BDRYX) 
is  shown  in  Fig.  8.  The  collision  density  and  tracklength  estimates 
of  the  flux  are  plotted  at  the  midpoint  between  two  boundaries.  The 
agreement  is  excellent  up  to  190  cm. 

The  comparison  of  the  flux  calculated  by  Monte  Carlo  using  the 
BDRYX  estimator  with  that  calculated  by  ANISN  is  shown  in  Table  7  and 
Fig.  9.  The  Monte  Carlo  results  are  the  points  and  the  ANISN  result 
is  the  smooth  curve.  The  agreement  between  the  two  calculations  is 
quite  good  and  within  the  statistical  variations  of  the  Monte  Carlo 
results.  For  distances  less  than  200  cm  the  statistical  errors  were 
less  than  3%.  In  the  vicinity  of  the  pressure  vessel  wall  the 
fractional  standard  deviation  is  more  like  10-20^.  For  a  more  precise 
comparison  in  this  region  further  biasing  techniques  to  reduce  the 
variance  would  be  helpful.  One  can  say,  however,  that  the  agreement 
between  ANISN  and  Monte  Carlo  near  the  pressure  vessel  wall  is  within 
20%. 


32 


neutrons /cm  -sec (E  >  l.OMeV) 


DISTANCE  FROM  CORE  CENTERLINE  (cm) 


Fig.  9  -  Fast  flux  calculated  by  MORSE  and  ANISN  using  slab  models 


A  comparison  of  the  neutron  energy  spectra  at  the  inside  of  the 
pressure  vessel  wall  is  shown  in  Fig.  10.  Again,  on  the  average,  the 
agreement  between  ANISN  and  Monte  Carlo  is  within  the  statistical 
uncertainty  of  the  Monte  Carlo  results. 

VIII.  MORSE  MONTE  CARLO  CALCULATION  IN  THE  CYLINDRICAL  MOOEL 

A.  Cylindrical  model  geometry  and  neutron  source  parameters. 

Since  the  major  reactor  components  are  more  nearly 
cylinders  than  slabs,  a  better  approximation  would  be  to  calculate  the 
neutron  flux  using  cylindrical  shells  instead  of  slabs.  In  the 
cylindrical  geometry  Monte  Carlo  calculation  also  the  height  of  each 
cylinder  was  taken  to  be  20,000  cm  and  the  outer  radius  of  each 
annulus  is  the  same  as  in  Table  1.  Since  the  reactor  core  and 
shielding  components  are  cylindrically  symmetric,  the  flux  will  be 
uniform  on  any  cylindrical  surface.  That  is,  the  flux  averaged  over 
the  surface  is  the  flux  at  any  point  on  the  surface.  Again  we  may 
replace  the  volume  source  with  a  line  source  and  use  MORSE  to 
calculate  the  average  flux  on  cylinders  of  radius  R. 
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FLUX  DENSITY  (n/cm  -sec-MeV) 


ENERGY (MeV) 


Fig.  10  -  Neutron  energy  spectra  at  the  inside  of  the  pressure 
vessel  wall  calculated  by  MORSE  and  ANISN  using  slab  models 
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For  the  cylindrical  calculation  we  used  the  same  power 
density  as  in  the  slab  case.  Since  the  volume  of  an  annulus  increases 
with  radius,  however,  the  distribution  of  relative  power  for  the  two 
cases  are  different. 


The  flux  estimator  used  in  the  cylindrical  model  was  the 
boundary  crossing  estimator,  suitably  modified  for  cylindrical 
symmetry.  The  boundary  crossing  estimator  in  the  cylindrical  case  is 
given  by 


*(R)  '  = 


U. 

IcosT)' 


where  R  is  the  core  radius, 
c 


(10) 


The  MORSE  Monte  Carlo  results  in  cylindrical  geometry  are 
shown  compared  to  ANISN  results  in  slab  geometry  in  Fig.  11.  Within 
185  cm  from  the  core  center  there  is  little  difference  outside 
statistics  between  the  calculations.  However,  beyond  135  cm  the  MORSE 
cylindrical  results  average  are  about  30%  lower  than  ANISN.  As  far  as 
the  radiation  safety  aspects  are  concerned  ANISN  is  somewhat 
conservative  in  predicting  a  higher  value  of  the  neutron  flux  in  the 
slab  model.  ANISN  calculations  in  cylindrical  geometry  were  not 
available. 
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neutrons  /cm  -  sec  (E  >1.0  MeV ) 


Fig.  11  -  Fast  flux  calculated  by  MORSE  using  a  cylindrical  model 
and  ANISN  using  a  slab  model 


FLUX  DENSITY  ( n/cm  -sec-MeV) 


Fig.  12  -  Neutron  energy  spectra  calculated  by  MORSE  using  a 
cylindrical  model  and  AN1SN  using  a  slab  model 
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FLUX  DENSITY  (n/cm  -sec-MeV) 


Fig.  13  -  Comparison  of  neutron  energy  spectra  calculated  by  MORSE 
using  cylindrical  and  slab  models 
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A  comparison  of  the  neutron  energy  spectra  at  the  inner 
surface  of  the  pressure  vessel  wall  for  ANISN  calculated  in  the  slab 
model  and  MORSE  calculated  in  the  cylindrical  model  is  shown  in  Fig. 
12.  Comparison  for  MORSE  in  the  slab  and  cylindrical  models  is  made 
in  Fig.  13.  In  both  comparisons  no  difference  in  neutron  spectra  are 
found  at  least  within  the  statistical  uncertainties  of  the  MORSE 
calculations. 

IX.  CONCLUSIONS 

The  fast  neutron  flux  in  the  reactor  core  midplane  calculated  in 
slab  geometry  with  the  codes  MORSE  and  ANISN  were  in  agreement  to 
within  the  statistical  error  of  the  MORSE  calculation.  The  neutron 
energy  spectra  at  the  inner  surface  of  the  pressure  vessel  wall  were 
also  in  agreement  in  this  situation.  MORSE  results  in  cylindrical 
geometry  show  the  effect  of  the  1/r  dependence  of  the  flux  outside  the 
core;  the  normalized  neutron  spectrum  at  the  pressure  vessel  agreed 
within  statistical  error  with  that  calculated  in  slab  geometry. 
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Appendix  I 


SUBROUTINE  SOURCE  ( IG,U,V,W,X,Y,Z, MATE, MED, AG, ISOUR.ITSTR, 
1NGPQT3,DDF, ISBIAS.NMTG) 

DIMENSION  P0Wf48),PSP0Wf43),RAD(48l 
COMMON  WTS(l) 

DATA  I N IT /O/ 

IF(INIT.EQ.l)  GO  TO  100 
READ  f 5 ,63)  MRVISO,BTHETA,BF 
R EAD C 5,63)  IZBIAS,EMFP,RCORE 
WRITE (6, 64)  MRVISO,BTHETA,BF 

64  FORMAT  ( 10X,  "MRVISO=" ,  15 , 10X,  "BTHETA="  ,E10. 5, 10X,  "BF=''  ,E10. 5) 
WR I TE f 6 , 65 )  IZBIAS,EMFP,RCORE 

65  FORMAT  f 10X," IZBIAS=" ,15, 10X, "EMFP=" ,E10. 5, 10X, "RCORE=" ,E10.5) 
IF ( IZBIAS.EQ. 1 )  S I GT = 1 ./EMFP 

IF  f I  ZB  I AS . NE . 1 )  S I GT =0 . 0 
63  FORMAT  f I5,5X,4E10.5) 

CTHETA=COS f  BTHETA*6 . 28/360 . ) 

FREQ=1 ./( 1 ,+BF) 

READI5, 10)  fRADfI),I=I,48) 

R EAD ( 5 , 10)  (POW( I ) , 1  =  1 ,48) 

10  FORMAT (5(E1 1 . 5,4X) ) 

C  CONSTRUCT  SPATIAL  GROUP  CDF 
PSPOW ( 1 )=POWf 1 ) 

DO  12  1=2,48 

12  PSPOW(I)=POW(I)+PSPOWfI-l) 

PSMAX=PSP0W(48) 

PSPOWf  l)=POWfl)*EXPf-(RCORE-RAD(l))*SIGT) 

DO  5  1=2,48 

5  PSPOWf I )=POW( I )*EXPf -fRCORE-RADfl ) ) *S I GT ) +PSPOW f I - 1 ) 
EPSMAX=PSPOWf  48) 

DO  6  1=1,48 

6  PSPOWf I )=PSPOWf I )/PSPOWf 48) 

INIT  =1 

100  CONTINUE 

C  CHOOSE  SPATIAL  GROUP 

P=FLTRNFfO) 

DO  7  1=1,48 

IF  (P.LE. PSPOWf I))  GO  TO  8 

7  CONTINUE 

8  Z=RADf I ) 

C  ASSIGN  WATE 

WATE=EXP ( +  ( RCOR  E-RAD ( I ) )*SIGT)*(EPSMAX/PSMAX)*WATE 
C  CHOOSE  ENERGY  GROUP 

IF ( I  SOUR)  15,15,60 
15  WATE=WATE*DDF 

IFfISBIAS)  20,20,25 
20  NWT=2*NMTG 

GO  TO  30 
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25  NWT=3*NMTG 
30  R=FITRNF (R) 

DO  35  1=1 ,NGPQT3 
IF  (R-  WTSf I+NWT) )  40,40,35 
35  CONTINUE 
40  IG=I 

IF  f I3BIAS1  60,60,45 
45  IF ( I- 1 )  60,50,55 

50  WATE=WATE*WT$ ( 2*NMTG+ 1 ) /WTS ( 3*NMTG+ 1 ) 

GO  TO  60 

55  WATE=WATE* f  WTS ( 2*NMTG+I ) -WTS  f  2*NMTG+I - 1 ) ) / 

1 ( WTS ( 3*NMTG+I ) -WTS ( 3*NMTG+I - 1 ) ) 

60  CONTINUE 

IF  ( MR V I  SO- 1 1  61,62,61 
62  CONTINUE 

C  SELECT  DIRECTION  COSINES 
2000  R1=FLTRNF(R) 

R2=FLTRNF (R) 

X1=2.0*R1-1 .0 
XSQ=X1*X1 
YSQ=R2*R2 
D=XSQ+YSO 

IF (D- 1 .01  2010,2010,2000 
2010  C0FI=(XSQ-YSQ1/D 
SIFI=2.*X1*R2/D 
R1=FLTRNF (R) 

R2=FLTRNF (R) 

IFfRl-FREQl  2030,2020,2020 
2020  W=CTHETA+( 1 .0-CTHETA)*R2 

AWATE =0 . 5* (1 . 0-CTHETA ) * ( 1 . 0+BF ) /BF 
GO  TO  2040 

2030  W=(l ,0+CTHETA)*R2-l .0 

AWATE =0 . 5* (1 . O+CTHETA ) * ( 1 . +BF ) 

2040  SITH=SQRT(1.0-W*W) 

U=SITH*COFI 

V=SITH*SIFI 

WATE=WATE*AWATE 

61  RETURN 
END 
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ApDendix  II 


SUBROUTINE  BDRYX 

COMMON  /POET/  ND,NNE,NE,NT,NA,NRESP,NEX,NEXND,NEND,NDNR,NTNR,NTNE, 

1  NANE,NTNDNR,NTNEND,NANEND,LOCRSP,LOCXD,LOCIB,LOCCO,LOCT,LOCUD, 

2  LOCSD,LOCQE,LOCQT,LOCQTE,LOCOAE,LMAX,EFIRST,EGTOP 

COMMON  /NUTRON/  NAME,NAMEX, IG, IGO,NMED,MEDOLD,NREG,U,V,W,UOLD,VOLD 
1  , WOLD, X,Y,Z,XOLD,YOLD,ZOLD,WATE,OLDWT,WTBC,BLZNT,BLZON, AGE, OLDAGE 
COMMON  BC( 1 ) 

DIMENSION  NC(  1 ) 

EQUIVALENCE  (BC( 1 ) ,NC( 1 ) ) 

Z2=0.999*Z 
Z22= 1 .OOI*Z 
IA  =  LOCXD  +  3*ND 
DO  10  1=1, ND 
IA  =  IA  +  1 

IF ( Z2-BC CIA)')  20,20,10 
10  CONTINUE 
GO  TO  100 

20  IF ( Z22-BC f I A ) )  30,30,40 

30  CONTINUE 
GO  TO  100 
40  COS=W- 1 . E- 10 
ABC  =  ABS  (COS) 

IF  (COS)  50,60,50 
50  IF  (ABC-1.0001)  70,60,60 
60  CALL  HELP  (4HBDRX,1 ,1 ,1 ,1) 

CALL  ERROR 

70  IF  (ABC-0.01)  80,90,90 
80  ABC  =  0.005 
90  CON  =  WATE/ABC 

CALL  FLUXST ( I , I G, CON, AGE, COS, 01 

C  *  *  SWITCH  =  0  -  STORE  IN  ALL  RELEVANT  ARRAYS  EXCEPT  UD 

INN  =  LOCXD  +  6*ND  +  I 

C  *  *THIS  STORE  IS  IN  THE  FIRST  OF  THE  NEXND  ARRAYS  SET  ASIDE  BY  SCORIN 
NC(INN)  =  NC(INN)  +  1 
100  CONTINUE' 

ZABS=ABS(Z) 

IF ( ZABS-0 . 00 1 )  1,2,2 

1  CALL  ALBDO ( IG , U , V , W, WATE , NMED , NREG 1 

2  CONTINUE 
RETURN 
END 
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Appendix  III 


SUBROUTINE  GTMED  (MOGEOM.MOXSEC) 
MDXSEC* ( 3*MOGEOM) /2- ( ( MOGEOM/2 )*2 ) 
IF(MOGEOH.EQ.IOOO)  MDXSEC-1000 
RETURN 
ENO 
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